rm(list = ls(all.names = TRUE))
gc()
set.seed(211917)
options(scipen=999)

packages <-c("tidyverse","foreign","stargazer","naniar","xtable","scales","girdExtra")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)
 
setwd("PUT YOUR WORKING DIRECTORY HERE")

df <- foreign::read.dta("./data.dta")

#---------------#
# Data cleaning #
#---------------#

df$respondent_christian <- 1 - df$respondent_muslim

## declare covariates
df$sample = 1
df$q6_years_married[df$q6_years_married == -9] <- 0
treats <- c("sample", "respondent_muslim", "respondent_christian")
vars <- c("respondent_male", "q1_age", "respondent_muslim", "q5_convert", "q6_married", "q6_years_married", "q7_has_job", "q8_has_income", "q8_income_last_month", "q3_kamba", "q3_kikuyu", "q3_luo", "q3_somali", "enumerator_muslim", "enumerator_male")

## helper function
get_stats <- function(var, treat, data) {
  c(mean(data[, var][data[, treat] ==1 ], na.rm = T), sd(data[, var][data[, treat] ==1 ], na.rm = T))
}

## apply function
btable <- lapply(treats, function(y) t(sapply(vars, function(x) {
  get_stats(x, y, data = df)
})))

## to df, add n, colnames
btable <- data.frame(do.call(cbind, btable))

## convert percentages
perc <- c("respondent_male","respondent_muslim","enumerator_muslim","enumerator_male","q3_kamba","q3_kikuyu","q3_luo","q3_somali","q5_convert","q6_married","q7_has_job","q8_has_income")
btable[rownames(btable) %in% perc, ] <- btable[rownames(btable) %in% perc, ] * 100

## add varnames
btable$variable <- rownames(btable)

## add n obs
n <- sapply(vars, function(x) sum(!is.na(df[, x])))
btable <- cbind(n, btable)

## change colnames
cn <- paste0(rep(treats, each = 2), rep(c("_mean", "_sd"), 2))
colnames(btable)[c(-1, -ncol(btable))] <- cn
btable <- data.frame(round(btable[, -which(colnames(btable) == "variable")], 1))
btable

library(xtable)
xtable(btable, digits=1)

